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Abstract 

A mathematical model for poro-visco-plastic com- 
paction and pressure solution in porous sediments 
has been formulated using the Voigt-type rheological 
constitutive relation as derived from experimental 
data. The governing equations reduce to a non- 
linear hyperbolic heat conduction equation in the 
case of slow deformation where permeability is 
relatively high and the pore fluid pressure is nearly 
hydrostatic, while travelling wave exists in the 
opposite limit where over-pressuring occurs and the 
pore fluid pressure is almost quasi-lithostatic. Full 
numerical simulation using a finite element method 
agree well with the approximate analytical solutions. 
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1 Introduction 

Many physical properties such permeability, viscos- 
ity. Young's modulus and thermal conductivity vary 
with porosity or fluid content in sedements and min- 
erals. The porosity in turn depends on the de- 
formation and compaction state and can be calcu- 
lated from the compaction curves resulting from the 
proper compaction modelling [Audet and Fowler, 
1992] . Furthmore, compaction is also related to over- 
pressuring, mineral deposition, hydrocarbon genera- 
tion and oil migration in reservoir. Thus the cor- 
rect modelling of compaction is both of scientific 
importance as well as industrial interest. However, 
due to the nonlinear feature in the compaction pro- 
cess and difficulty in formulating the accurate and 
realistic rheological relationship in sediments and 
rocks, most of existing models use simplified rheol- 
ogy such as poroelastic, purely viscous or viscoelastic 
relations [Rutter, 1976; Wangen, 1992; Holzbercher, 
1998; Revil, 1999; Yang,2000]. The rheological prop- 
erties of realistic granular sediments are usually vis- 
coelastic or viscoplastic as implied by experiments. 



Yang [2000] presents a viscoelastic model of Maxwell 
type and comparison of analytical solutions with the 
numerical simulations shows very good agreement. 
However, ReviTs [1999] work suggests that it maybe 
more appropriate to use a poro-visco-plastic model 
of Voigt type. We intend to do the analysis similar 
to Yang's [2000] but using the Voigt-type consitutive 
relation as derived by Revil [1999]. 

Much of the work in this area has been re- 
viewed by Rieke and Chilingarian [1974], Birch- 
wood and Turcotte [1994] and later by Fowler and 
Yang[1999]. This paper aims at providing a new ap- 
proach to compaction and pressure solution by us- 
ing ReviTs visco-poro-plastic relation of Voigt type 
[Revil, 1999]. The nonlinear partial differential equa- 
tions are then analysed by using asymptotic methods 
and the obtained analytical solutions are compared 
with numerical simulations. Although the present 
work mainly concerns the 1-D theorectical formula- 
tion and analytical solution procedure, however, we 
intend to provide a simplified and yet realistic frame- 
work for further research in this area and shows how 
compaction mechanism is related to rheological re- 
lationships and material properties of porous sedi- 
ments, so that more realistic constitutive relation- 
ships can be formulated and analysed. 

2 MATHEMATICAL MODEL 

The fundamental model of compaction and pressure 
solution is essentially similar to the model of soil con- 
solidation process. The solid sediments act as a com- 
pressible porous matrix, so mass conservation of pore 
fiuid together with Darcy's law leads to an equation 
of the general type. Based on earlier work by Au- 
det and Fowler [1992], Revil[im9\ and Yanc/[2000], 
we can write down the poro-visco-plastic compaction 
model of Voigt type. In the one-dimensional case, we 
have the following governing equations 
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where (f) is porosity, and are the velocities of 
fluid and solid matrix, respectively, k and /i are the 
matrix permeability and the liquid viscosity, pi and 
Ps are the densities of fluid and solid matrix, is the 
effective pressure, is the pore pressure, and g is the 
gravitational acceleration. G = 1 +4?7o/3^o is a con- 
stant describing the material properties with r}o and 
^0 being the shear modulus and bulk viscosity [Bird 
et al, 1977]. The first two equations are the conser- 
vation of mass for the solid phase and liquid phase, 
respectively. The third equation is the Darcy's law in 
1-D form and the last equation is actually the force 
balance in a simplified form whose detailed deriva- 
tion can be found in [Fowler and Yang, 1998]. Com- 
bining equation (3) and (4), we have 
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In writing the above equations, we have used an up- 
ward coordinate z originating from z = 0, which 
corresponds to the bottom of the sedimentary col- 
umn, so that the ocean fioor z = h{t) moves as com- 
paction proceeds. We use such a coordinate system 
because it simplifies the analytical solution proce- 
dure and also in keeping with the similar lines of 
earlier work in this area [Audet and Fowler, 1992; 
Yang, 2000]. However, the conventional depth co- 
ordinate is simply z — h{t), thus the transformation 
shall be straightforward once the basin thickness h{t) 
is known. As we shall see in the later sections, we 
provide an explicit formula for h{t) as a very good 
approximation. 

In addition, a rheological compactional relation- 
ship derived from experimental data [Revil, 1999] is 
needed to complete this model in the form 



Pe = -<e^ - E / —dt, 

dz Jo dz 



(6) 



where E is the elastic modulus and ^ is the viscosity 
modulus. There are essentially the same parame- 
ters as introduced by Revil [1999]. The first term of 
the right hand of the equation is the usual contribu- 
tion by viscous plastic deformation, while the second 
term corresponds to the poro-elastic deformation. 

3 Non-dimensionalization 

To write the governing equations in dimensionless 
forms, typical length and time scales are required. 



For a typical sedimentation rate rhs, the correspond- 
ing typical time scale is rf/m,, where the typical 
length scale d can be defined as 
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so that the dimensionless pressure p = Gpe/{ps — 
Pi)gd = 0{1). Meanwhile, we scale z with d, 
with rhs, time t with d/rfis, permeability k with kg. 
By writing k{(j)) = kok* , z = dz*, and dropping 
the asterisks, we thus have 
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Adding (8) and (9) together and integrating from 
the bottom, we have 
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where u = (f>{u'' — u^) is the Darcy flow velocity. Now, 
we have 

u^ = x{^r[-^-{i -<!>)]. (15) 

00 dz 
du' „ /■* du" , 

where we have used the nonlinear constitutive re- 
lation for permeability k{(j)) of typical form [Smith, 
1971] 

m = i^r, (17) 

where 0o is the initial depositional porosity. The 
exponent m has a typical value of 3 6 for sands 
and sandstones. 

The boundary conditions are 



■0) = (or equivalently, 



p = 0. 



0), 



at z = 0, 
(18) 



2 



nis 



■A(-fr 



dp 



(1 -(!>)] at z^h{t). (19) 



It is useful to estimate these parameters by using 
values taken from observations. By using the typi- 
cal values of pi ~ lO^kgm"^, ps ~ 2.5 x lO^kgm"^, 
fco - 10-1^-10-2" m2, /i ~ 10-3Nsm^ ^ - IxlO^i 
N s m-2, m, ~ 300 m Ma-^ = 1 x 10-" m s"\ g - 
lOms-2, E ~ lO^N/m^ G - 1, d ~ 1000m ; then 
A « 0.01 ~ 1000 and 5 40. Wo can sec that the 
main parameters A and 5, which govern the evolu- 
tion of the fluid flow and porosity in sedimentary 
basins, are the ratios of permeability to sedimenta- 
tion rate and material mochilus to the typical pres- 
sure scale. As A is essentially controlled by the hy- 
draulic conductivity, so it becomes the dominant pa- 
rameter controlling the whole compaction and pres- 
sure solution processes. 



4 Asymptotic Analysis 

Since the nondimensional parameter A ~ 0.01 ~ 
1000 varies greatly and essentially controls the com- 
paction process, we can expect that the two distin- 
guished limits (A ^ 1 and A » 1) will have very 
different features in porosity and flow evolutions. In 
fact, A = 1 dcflnes a transition between slow com- 
paction (A <C 1) and fast compaction (A ^ 1). The 
case of A <C 1 corresponds to the situation where 
the pore fluid pressure is nearly hydrostatic whereas 
the opposite case corresponds to an overpressured 
section in which the pore fluid pressure is quasi- 
lithostatic. Thus we can follow the similar asymp- 
totic analysis [Fowler and Yang, 1998,1999] to obtain 
some analytical asymptotic solutions. 



combining these above three equations, we have a 
single equation for (j) 
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As the s ^ 1 and A ^ 1, we can use the approxi- 
mation 04 « AS(/)22 so that we have 
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This problem is in fact equivalent to the problem 

of hyperbolic heat condiiction or non-Fourier heat 
equation which is well-documented in heat transfer 
and laser pulse modelling [Antaki, 1997]. By using 
the Laplace transform method, wc can write the so- 
lution approximately in terms of Bessel functions Jj 

as 
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and ai is the ith real non-negative root of equation 
Ji{cti) = 0. We can see that compaction essentially 
occurs in a boundary layer near the bottom with a 
thickness of the order of VX. 



4.1 Slow Deformation (A <C 1) 

In the nearly hydrostatic case ofA^l,z^l,f^l, 
p 1 implies that <C 1 and §f ~ 0, then (p « (pg. 
We thus have 
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4.2 Fast Deformation (A > 1) 

In the case of A 1, the dependence of permeability 
on porosity {(p/(f)o)"^ decrease dramatically, so that 
A(0/(/)o)™ is only bigger enough when > 0* = 
(/>oexp[— (InA)/™]. Thus, we have 
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and using the equation (29), we get 
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Combining these equation, we can get a single equa- 
tion for 6 
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Now wc can seek the travehng wave solution of the 
form <f) = (f){C,) with Q = z — ct, so that we have 



(1 -0)(1 = -c(/." + S/' 



(33) 



where 0' = dcji/dC,. We can easily write the solution 
as 



*^l-(l-«exp[H^v5^<i^|(34) 

In fact, the above solution is only valid for the top 
part when (f> < (f>^. Using equation (6) and (j) ^ (j)^ 
as C — >■ —00, the travelling wave implies that 

c0 + (1 - (j))u' = C(/>o + (to, - c)(l - (j)o). (35) 

so that we have 



c i=a m. 



which means the basin thickness increases linearly 
with time. 



5 Numerical Simulations 

In order to check the accuracy of the above analysis, 
we used a finite element method to solve the above 
equations (8)-(ll). For simplicity, we only present 
the related results in Fig.l where t = 5 with different 
values of the A = 0.001,0.1,10,1000, 0o = 0.5 and 
S = 40. This good agreement near the top and the 
bottom regions suggest that there exists a travelling 
wave solution on the top for A » 1 and a boundary 
layer near the bottom for the small A case. 



the bottom of the compacting column. On the other 

hand, for the large A ^ 1 case where sedimentation 
or loading is slow in high permeable sediments, the 
travelling wave solution exists and the top surface 
moves up with nearly constant velocity. 

Compared with the earlier work [Yang, 2000] 
using the viscoelastic rheological relation, we see 
that the there is no essential difference in the case 
small deformation (A <C 1). Boundary layer exists 
in both viscoelastic and poro-viscous-plastic cases 
although the slight difference is that the former 
viscoelastic case corresponds to the heat conduction 
with a constant flux and a constant source term, 
while the latter poro-visco-plastic case is mainly a 
hyperbolic heat condution mechanism. The case 
of fast deformation and compaction (A » 1) is 
more complicated. Although both viscoelastic and 
poro-viscous-plastic cases have a transition at the 
depth where ^ « however, the mechanisms 
above and below the transition are very different. 
In the viscoelastic case, the top region is essentially 
poroelastic while the lower region is almost purely 
viscous, while in the present poro-visco-plastic case, 
the mechanism is a complicated combination of 
viscous-plastic mechanism and poro-elastic defor- 
mation process controlled by the proper balance 
of the parameters A and S. Furthermore, in com- 
parison with the earlier results, this work suggests 
that Voigt-type poro-visco-plastic deformation 
has much interesting characteristics due to its 
time-dependence feature. 
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6 Discussion 

A mathematical model for poro-visco-plastic 
comapction and pressure solution in porous sed- 
iments has been formulated using the Voigt-type 
rheological constitutive relation as derived from 
experimental data by Revil [1999]. After the proper 
scalings, the governing equations reduce to a system 
of coupled partial diff'erential equations of mixed 
type. In the case of small A where the sedimentation 
is fast, permeability is small and the pore fluid 
pressure is nearly hydrostatic, the pressure solution 
process reduces to the case of hyperbolic heat con- 
duction equation with a boundary layer forming at 
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